This code reproduces the trial-level analyses reported in the following manuscript:

Cosme, D., Mobasser, A., & J. H. Pfeifer. If you’re happy and you know it: Neural correlates of self-evaluated psychological health and well-being

load packages

if(!require('pacman')) {
    install.packages('pacman')
}

pacman::p_load(tidyverse, gtools, GGally, sjstats, lme4, lmerTest, knitr, ggeffects, kableExtra, psych, install = TRUE)

define aesthetics

rois = c("#006989", "#56445D", "#8EC922")
constructs = c("#FEC601", "#254E70", "#F43C13")
instructions = wesanderson::wes_palette("Darjeeling1", 2, "continuous")
valence = c("#119DA4", "#19647E")
plot_aes = theme_minimal() +
  theme(legend.position = "top",
        legend.text = element_text(size = 12),
        text = element_text(size = 16, family = "Futura Medium"),
        axis.text = element_text(color = "black"),
        axis.line = element_line(colour = "black"),
        axis.ticks.y = element_blank(),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        panel.border = element_blank(),
        panel.background = element_blank())

define functions

model_table = function(model) {
  model %>%
    broom.mixed::tidy(., conf.int = TRUE) %>%
    filter(effect == "fixed") %>%
    rename("SE" = std.error,
           "z" = statistic,
           "p" = p.value) %>%
    mutate(term = gsub("\\(Intercept\\)", "Intercept (self-oriented well-being)", term),
           term = gsub("log\\(RT\\)", "RT", term),
           term = gsub("constructill-being", "Construct (ill-being)", term),
           term = gsub("constructsocial well-being", "Construct (social well-being)", term),
           term = gsub(":", " x ", term),
           p = ifelse(p < .001, "< .001",
                      ifelse(p == 1, "1.000", gsub("0.(.*)", ".\\1", sprintf("%.3f", p)))),
           `b [95% CI]` = sprintf("%.2f [%0.2f, %.2f]", estimate, conf.low, conf.high),
           z = abs(round(z, 2))) %>%
    select(term, `b [95% CI]`, z, p) %>%
    kable() %>%
    kableExtra::kable_styling()
}

load data

task = read.csv("../../data/task_data.csv") %>%
  mutate(responseYN = as.factor(responseYN),
         construct = recode(construct, "well-being" = "self-oriented well-being",
                            "social" = "social well-being"),
         construct = factor(construct, levels = c("self-oriented well-being", "social well-being", "ill-being")),
         valence = factor(valence, levels = c("positive", "negative")))
betas = read.csv("../../data/neuro_data.csv")

exclude outliers and standardize

  • Subset self trials
  • Exclude outlier trials that are > 3 SDs from roi median across participants
trial_conditions = task %>%
  select(subjectID, trial, run, instruction)

betas_outlier = betas %>%
  left_join(., trial_conditions) %>%
  filter(instruction == "self") %>%
  group_by(roi) %>%
  mutate(median = median(meanPE, na.rm = TRUE),
         sd3 = 3*sd(meanPE, na.rm = TRUE),
         outlier = ifelse(meanPE > median + sd3 | meanPE < median - sd3, "yes", "no")) 
## Joining, by = c("subjectID", "trial", "run")
betas_outlier %>%
  group_by(outlier) %>%
  summarize(n = n()) %>%
  ungroup() %>%
  mutate(total = sum(n),
         percent = round((n / total) * 100, 2)) %>%
  kable(format = "pandoc")
outlier n total percent
no 11900 11988 99.27
yes 88 11988 0.73
betas_ex = betas_outlier %>%
  filter(outlier == "no") %>%
  select(-c(median, sd3, outlier))

merge data and exclude participants

  • Motion exclusions: FP091
  • Technical failure: FP080, FP082
  • Non-compliance: FP021, FP049, FP085, FP121
  • Standardize within participant and ROI
data_ex = task %>%
  left_join(., betas_ex) %>%
  filter(!subjectID %in% c("FP091", "FP080", "FP082", "FP021", "FP049", "FP085", "FP121")) %>%
  filter(!is.na(meanPE)) %>%
  filter(instruction == "self") %>%
  group_by(roi, subjectID) %>%
  mutate(std = meanPE / sd(meanPE, na.rm = TRUE)) %>%
  select(-sdPE)
## Joining, by = c("trial", "subjectID", "run", "instruction")

remove missing data for MLM analysis

data_complete_mod = data_ex %>%
  select(-meanPE) %>%
  spread(roi, std) %>%
  select(subjectID, trial, run, construct, item, valence, response, responseYN, RT_original, RT, pgACC, vmPFC, VS) %>%
  na.omit()

item responses

percent endorsements

percents = data_complete_mod %>%
  group_by(construct, item, response) %>%
  summarize(n = n()) %>%
  spread(response, n) %>%
  mutate(no = ifelse(is.na(no), 0, no),
         yes = ifelse(is.na(yes), 0, yes),
         total = no + yes) %>%
  gather(response, n, no, yes) %>%
  mutate(percent = round((n/total) * 100, 1),
         n = sprintf("%s / %s", n, total)) %>%
  filter(response == "yes")
## `summarise()` has grouped output by 'construct', 'item'. You can override using
## the `.groups` argument.
percents %>%
  select(-response, -total) %>%
  arrange(desc(percent))  %>%
  kable() %>%
  kableExtra::kable_styling()
construct item n percent
self-oriented well-being capable 5 / 5 100.0
social well-being have support systems 5 / 5 100.0
social well-being well-liked 5 / 5 100.0
ill-being tired 5 / 5 100.0
self-oriented well-being interested in life 97 / 98 99.0
social well-being loved 102 / 103 99.0
social well-being trusted 101 / 104 97.1
social well-being am cared for 98 / 102 96.1
social well-being comfortable with friends 96 / 104 92.3
self-oriented well-being happy 94 / 102 92.2
self-oriented well-being live a meaningful life 95 / 103 92.2
self-oriented well-being achieve goals 95 / 104 91.3
social well-being have warm relationships 95 / 104 91.3
self-oriented well-being hopeful 88 / 97 90.7
self-oriented well-being live a purposeful life 92 / 102 90.2
self-oriented well-being engaged in tasks 87 / 97 89.7
social well-being have good friends 90 / 102 88.2
social well-being have trusting relationships 87 / 99 87.9
self-oriented well-being joyful 88 / 103 85.4
self-oriented well-being like myself 84 / 103 81.6
social well-being feel encouraged 82 / 101 81.2
self-oriented well-being positive 79 / 98 80.6
self-oriented well-being full of energy 4 / 5 80.0
self-oriented well-being optimistic 4 / 5 80.0
social well-being helped by friends 76 / 96 79.2
self-oriented well-being satisfied with life 79 / 101 78.2
self-oriented well-being proud 78 / 100 78.0
ill-being stressed 67 / 100 67.0
social well-being belong in social group 62 / 95 65.3
self-oriented well-being cheerful 3 / 5 60.0
ill-being anxious 57 / 104 54.8
ill-being nervous 53 / 103 51.5
ill-being overwhelmed 49 / 103 47.6
social well-being lonely 31 / 102 30.4
ill-being fearful 26 / 98 26.5
ill-being unhappy 1 / 4 25.0
social well-being isolated 24 / 102 23.5
ill-being feel unimportant 18 / 103 17.5
ill-being feel rejected 17 / 104 16.3
ill-being sad 14 / 99 14.1
ill-being angry 10 / 98 10.2
ill-being feel useless 9 / 99 9.1
ill-being unable to cope 8 / 104 7.7
ill-being feel worthless 5 / 99 5.1
social well-being alienated 0 / 5 0.0
social well-being neglected 0 / 5 0.0
ill-being depressed 0 / 4 0.0
ill-being feel invisible 0 / 5 0.0

percent endorsements for improperly randomized items

percents %>%
  filter(total < 10) %>%
  select(-response, -total) %>%
  arrange(desc(construct)) %>%
  kable() %>%
  kableExtra::kable_styling()
construct item n percent
ill-being depressed 0 / 4 0
ill-being feel invisible 0 / 5 0
ill-being tired 5 / 5 100
ill-being unhappy 1 / 4 25
social well-being alienated 0 / 5 0
social well-being have support systems 5 / 5 100
social well-being neglected 0 / 5 0
social well-being well-liked 5 / 5 100
self-oriented well-being capable 5 / 5 100
self-oriented well-being cheerful 3 / 5 60
self-oriented well-being full of energy 4 / 5 80
self-oriented well-being optimistic 4 / 5 80

plot

data_complete_mod %>%
  group_by(construct, item, responseYN) %>%
  summarize(n = n()) %>%
  spread(responseYN, n) %>%
  mutate(no = ifelse(is.na(no), 0, no),
         yes = ifelse(is.na(yes), 0, yes),
         total = no + yes) %>%
  gather(responseYN, n, no, yes) %>%
  mutate(percent = round((n/total) * 100, 1),
         n = sprintf("%s / %s", n, total)) %>%
  filter(responseYN == "yes") %>%
  ggplot(aes(percent, fill = construct)) +
  geom_density(color = NA, alpha = .8) +
  scale_fill_manual(values = constructs) +
  labs(x = "\n percent endorsed") +
  plot_aes
## `summarise()` has grouped output by 'construct', 'item'. You can override using
## the `.groups` argument.

well-being facet correlations

data_corr = data_complete_mod %>%
  mutate(responseNum = recode(responseYN, "yes" = 1, "no" = 0),
         construct = recode(construct, "self-oriented well-being" = "self",
                            "social well-being" = "social")) %>%
  group_by(construct, subjectID) %>%
  summarize(sum = sum(responseNum, na.rm = TRUE),
            n = n()) %>%
  mutate(percent = (sum / n) * 100) %>%
  select(subjectID, construct, percent) %>%
  spread(construct, percent)
## `summarise()` has grouped output by 'construct'. You can override using the
## `.groups` argument.
# correlations
cor.ci(data_corr[,2:4])

## Call:corCi(x = x, keys = keys, n.iter = n.iter, p = p, overlap = overlap, 
##     poly = poly, method = method, plot = plot, minlength = minlength, 
##     n = n)
## 
##  Coefficients and bootstrapped confidence intervals 
##           self  socil ill-b
## self       1.00            
## social     0.57  1.00      
## ill-being -0.39 -0.42  1.00
## 
##  scale correlations and bootstrapped confidence intervals 
##             lower.emp lower.norm estimate upper.norm upper.emp p
## self-socil       0.42       0.42     0.57       0.71      0.71 0
## self-ill-b      -0.50      -0.53    -0.39      -0.24     -0.26 0
## socil-ill-b     -0.60      -0.59    -0.42      -0.22     -0.25 0

ROI descriptives

means, SDs, and correlations

# correlations
pgACC_vmPFC = data_complete_mod %>%
   mutate(subjectID = as.factor(subjectID)) %>%
   select(subjectID, pgACC, vmPFC) %>%
   rmcorr::rmcorr(subjectID, pgACC, vmPFC, .)

pgACC_VS = data_complete_mod %>%
   mutate(subjectID = as.factor(subjectID)) %>%
   select(subjectID, pgACC, VS) %>%
   rmcorr::rmcorr(subjectID, pgACC, VS, .)

vmPFC_VS = data_complete_mod %>%
   mutate(subjectID = as.factor(subjectID)) %>%
   select(subjectID, vmPFC, VS) %>%
   rmcorr::rmcorr(subjectID, vmPFC, VS, .)

# means and SDs
mean_table = data_complete_mod %>%
  gather(roi, value, pgACC, vmPFC, VS) %>%
  group_by(roi) %>%
  summarize(M = round(mean(value, na.rm = TRUE), 2),
            SD = round(sd(value, na.rm = TRUE), 2)) %>%
  mutate(order = ifelse(roi == "pgACC", 1,
                 ifelse(roi == "vmPFC", 2, 3)),
         roi = ifelse(roi == "pgACC", "1. pgACC",
               ifelse(roi == "vmPFC", "2. vmPFC", "3. VS"))) %>%
  rename("ROI" = roi) %>%
  arrange(order) %>%
  select(-order)

# table
corr_table = data.frame(ROI = c("1. pgACC", "2. vmPFC", "3. VS"),
           `1` = c("--", paste0(round(pgACC_vmPFC[[1]],2), " [",round(pgACC_vmPFC[[4]][1],2),", ",round(pgACC_vmPFC[[4]][2],2),"]"),
                   paste0(round(pgACC_VS[[1]],2), " [",round(pgACC_VS[[4]][1],2),", ",round(pgACC_VS[[4]][2],2),"]")),
           `2` = c("", "--", paste0(round(vmPFC_VS[[1]],2), " [",round(vmPFC_VS[[4]][1],2),", ",round(vmPFC_VS[[4]][2],2),"]")),
           `3` = c("", "", "--"),
           check.names = FALSE)

mean_table %>%
  full_join(., corr_table) %>%
  knitr::kable(format = "pandoc", caption = "Repeated Measures Correlations Between ROIs (n = 3694). All correlations are statistically significant, p < .001. 95% confidence intervals are bracketed")
## Joining, by = "ROI"
Repeated Measures Correlations Between ROIs (n = 3694). All correlations are statistically significant, p < .001. 95% confidence intervals are bracketed
ROI M SD 1 2 3
1. pgACC 0.04 1.03 –
2. vmPFC 0.15 1.02 0.84 [0.83, 0.85] –
3. VS 0.21 1.03 0.42 [0.39, 0.45] 0.45 [0.43, 0.48] –

RT descriptives

across facets and responses

data_complete_mod %>%
  ungroup() %>%
  summarize(`mean RT` = mean(RT_original, na.rm = TRUE),
            `median RT` = median(RT_original, na.rm = TRUE),
            SD = sd(RT_original, na.rm = TRUE)) %>%
    kable(digits = 2) %>%
    kableExtra::kable_styling()
mean RT median RT SD
1.62 1.48 0.61

across responses

data_complete_mod %>%
  group_by(construct) %>%
  summarize(`mean RT` = mean(RT_original, na.rm = TRUE),
            `median RT` = median(RT_original, na.rm = TRUE),
            SD = sd(RT_original, na.rm = TRUE)) %>%
    kable(digits = 2) %>%
    kableExtra::kable_styling()
construct mean RT median RT SD
self-oriented well-being 1.53 1.40 0.57
social well-being 1.61 1.44 0.62
ill-being 1.72 1.58 0.62
data_complete_mod %>%
  ggplot(aes(RT_original, fill = construct)) +
  geom_density(color = NA, alpha = .75) +
  scale_fill_manual(values = constructs) +
  labs(x = "\nRT (s)") +
  plot_aes

by response

data_complete_mod %>%
  group_by(construct, responseYN) %>%
  summarize(`mean RT` = mean(RT_original, na.rm = TRUE)) %>%
  spread(responseYN, `mean RT`) %>%
    kable() %>%
    kableExtra::kable_styling()
## `summarise()` has grouped output by 'construct'. You can override using the
## `.groups` argument.
construct no yes
self-oriented well-being 2.015276 1.455663
social well-being 2.113079 1.530455
ill-being 1.709059 1.746566

model testing diffferences

model_rt = lmer(RT_original ~ responseYN * construct + (1 | subjectID),
                data = data_complete_mod)

summary(model_rt)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: RT_original ~ responseYN * construct + (1 | subjectID)
##    Data: data_complete_mod
## 
## REML criterion at convergence: 5910.1
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.8363 -0.6517 -0.1830  0.4757  5.8536 
## 
## Random effects:
##  Groups    Name        Variance Std.Dev.
##  subjectID (Intercept) 0.07652  0.2766  
##  Residual              0.26907  0.5187  
## Number of obs: 3694, groups:  subjectID, 105
## 
## Fixed effects:
##                                            Estimate Std. Error         df
## (Intercept)                                 1.93039    0.05072  894.26616
## responseYNyes                              -0.46016    0.04614 3632.65369
## constructsocial well-being                  0.08782    0.05777 3599.67437
## constructill-being                         -0.18627    0.04660 3629.89748
## responseYNyes:constructsocial well-being   -0.00823    0.06222 3602.24250
## responseYNyes:constructill-being            0.38339    0.05826 3646.40501
##                                          t value             Pr(>|t|)    
## (Intercept)                               38.062 < 0.0000000000000002 ***
## responseYNyes                             -9.973 < 0.0000000000000002 ***
## constructsocial well-being                 1.520                0.129    
## constructill-being                        -3.997      0.0000653894996 ***
## responseYNyes:constructsocial well-being  -0.132                0.895    
## responseYNyes:constructill-being           6.581      0.0000000000535 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) rspnYN cnwll- cnstr- rYN:w-
## responsYNys -0.794                            
## cnstrctwll- -0.607  0.670                     
## cnstrctll-b -0.784  0.870  0.661              
## rspnsYNy:w-  0.566 -0.712 -0.932 -0.616       
## rspnsYNys:-  0.643 -0.810 -0.532 -0.823  0.565
predicted_sub = ggpredict(model_rt, terms = c("responseYN", "construct", "subjectID"), type = "random") %>%
  data.frame()

ggeffect(model_rt, terms = c("responseYN", "construct")) %>%
  data.frame() %>%
  ggplot(aes(x, predicted)) +
  #geom_line(data = predicted_sub, aes(color = group, group = interaction(facet, group)), size = .25, alpha = .25) +
  geom_line(aes(color = group, group = group), size = 1.5) +
  geom_pointrange(aes(color = group, ymin = conf.low, ymax = conf.high)) +
  scale_color_manual(values = constructs, name = "") + 
  scale_fill_manual(values = constructs, name = "") + 
  labs(x = "\nresponse", y = "reaction time (s)\n") + 
  plot_aes
## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

RT models

These are the primary analyses reported in the main text of the manuscript, controlling for reaction time.

run models

compare models

Determine the best fitting model.

model_0_rt = Base model including including no ROIs

model_1_rt = Main effects model including ROIs

model_2_rt = Interaction model including ROIs and their interactions with well-being construct

model_3_rt = Interaction model including ROIs and their interactions with valence rather than well-being construct

model_0_rt = glmer(responseYN ~ construct * RT + 
                     (1 | subjectID), 
                family = binomial, 
                data = data_complete_mod, control = glmerControl(optimizer = "bobyqa"))

model_1_rt = glmer(responseYN ~ construct * RT + pgACC + vmPFC + VS + 
                     (1 | subjectID), 
                family = binomial, 
                data = data_complete_mod, control = glmerControl(optimizer = "bobyqa"))

model_2_rt = glmer(responseYN ~ RT*construct + pgACC*construct + vmPFC*construct + VS*construct +
                     (1 | subjectID), 
                family = binomial, 
                data = data_complete_mod, control = glmerControl(optimizer = "bobyqa"))

model_3_rt = glmer(responseYN ~ RT*valence + pgACC*valence + vmPFC*valence + VS*valence +
                     (1 | subjectID), 
                family = binomial, 
                data = data_complete_mod, control = glmerControl(optimizer = "bobyqa"))

# compare models
anova(model_0_rt, model_1_rt, model_2_rt, model_3_rt)

calculate R2

Compare the model R2 for the best fitting compared to the base model.

performance::r2(model_0_rt)
## # R2 for Mixed Models
## 
##   Conditional R2: 0.512
##      Marginal R2: 0.479
performance::r2(model_2_rt)
## # R2 for Mixed Models
## 
##   Conditional R2: 0.532
##      Marginal R2: 0.498

check variance inflation in the best fitting model

# vif
car::vif(model_2_rt)
##                      GVIF Df GVIF^(1/(2*Df))
## RT               3.985384  1        1.996343
## construct        1.514415  2        1.109331
## pgACC           13.575127  1        3.684444
## vmPFC           13.404569  1        3.661225
## VS               4.932306  1        2.220880
## RT:construct     4.236780  2        1.434693
## construct:pgACC 40.780615  2        2.527048
## construct:vmPFC 48.147820  2        2.634172
## construct:VS     6.643231  2        1.605443

model summary

summary(model_2_rt)
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: binomial  ( logit )
## Formula: responseYN ~ RT * construct + pgACC * construct + vmPFC * construct +  
##     VS * construct + (1 | subjectID)
##    Data: data_complete_mod
## Control: glmerControl(optimizer = "bobyqa")
## 
##      AIC      BIC   logLik deviance df.resid 
##   3082.3   3181.7  -1525.1   3050.3     3678 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -8.1633 -0.5202  0.2040  0.3958  3.1031 
## 
## Random effects:
##  Groups    Name        Variance Std.Dev.
##  subjectID (Intercept) 0.2382   0.488   
## Number of obs: 3694, groups:  subjectID, 105
## 
## Fixed effects:
##                                  Estimate Std. Error z value
## (Intercept)                       1.89473    0.11189  16.933
## RT                               -3.15503    0.29977 -10.525
## constructsocial well-being        0.08607    0.14486   0.594
## constructill-being               -2.95012    0.12509 -23.584
## pgACC                             0.44782    0.18871   2.373
## vmPFC                            -0.39372    0.18837  -2.090
## VS                                0.10966    0.11200   0.979
## RT:constructsocial well-being    -0.03285    0.39988  -0.082
## RT:constructill-being             3.16002    0.35531   8.894
## constructsocial well-being:pgACC  0.06793    0.25779   0.264
## constructill-being:pgACC         -0.38844    0.22428  -1.732
## constructsocial well-being:vmPFC -0.23103    0.26387  -0.876
## constructill-being:vmPFC          0.75862    0.22610   3.355
## constructsocial well-being:VS    -0.03716    0.15559  -0.239
## constructill-being:VS            -0.18072    0.13704  -1.319
##                                              Pr(>|z|)    
## (Intercept)                      < 0.0000000000000002 ***
## RT                               < 0.0000000000000002 ***
## constructsocial well-being                   0.552418    
## constructill-being               < 0.0000000000000002 ***
## pgACC                                        0.017643 *  
## vmPFC                                        0.036604 *  
## VS                                           0.327535    
## RT:constructsocial well-being                0.934521    
## RT:constructill-being            < 0.0000000000000002 ***
## constructsocial well-being:pgACC             0.792143    
## constructill-being:pgACC                     0.083287 .  
## constructsocial well-being:vmPFC             0.381276    
## constructill-being:vmPFC                     0.000793 ***
## constructsocial well-being:VS                0.811232    
## constructill-being:VS                        0.187246    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation matrix not shown by default, as p = 15 > 12.
## Use print(x, correlation=TRUE)  or
##     vcov(x)        if you need it

model table

model_table(model_2_rt)
term b [95% CI] z p
Intercept (self-oriented well-being) 1.89 [1.68, 2.11] 16.93 < .001
RT -3.16 [-3.74, -2.57] 10.52 < .001
Construct (social well-being) 0.09 [-0.20, 0.37] 0.59 .552
Construct (ill-being) -2.95 [-3.20, -2.70] 23.58 < .001
pgACC 0.45 [0.08, 0.82] 2.37 .018
vmPFC -0.39 [-0.76, -0.02] 2.09 .037
VS 0.11 [-0.11, 0.33] 0.98 .328
RT x Construct (social well-being) -0.03 [-0.82, 0.75] 0.08 .935
RT x Construct (ill-being) 3.16 [2.46, 3.86] 8.89 < .001
Construct (social well-being) x pgACC 0.07 [-0.44, 0.57] 0.26 .792
Construct (ill-being) x pgACC -0.39 [-0.83, 0.05] 1.73 .083
Construct (social well-being) x vmPFC -0.23 [-0.75, 0.29] 0.88 .381
Construct (ill-being) x vmPFC 0.76 [0.32, 1.20] 3.36 < .001
Construct (social well-being) x VS -0.04 [-0.34, 0.27] 0.24 .811
Construct (ill-being) x VS -0.18 [-0.45, 0.09] 1.32 .187

marginal probabilities

all

ggeffect(model_2_rt)
## $RT
## # Predicted probabilities of responseYN
## 
##    RT | Predicted |       95% CI
## --------------------------------
## -1.50 |      0.98 | [0.97, 0.99]
## -1.00 |      0.95 | [0.94, 0.97]
## -0.50 |      0.88 | [0.85, 0.90]
##  0.00 |      0.72 | [0.69, 0.74]
##  0.50 |      0.47 | [0.42, 0.52]
##  1.00 |      0.23 | [0.18, 0.30]
## 
## $construct
## # Predicted probabilities of responseYN
## 
## construct                | Predicted |       95% CI
## ---------------------------------------------------
## self-oriented well-being |      0.90 | [0.88, 0.92]
## social well-being        |      0.91 | [0.88, 0.93]
## ill-being                |      0.27 | [0.24, 0.30]
## 
## $pgACC
## # Predicted probabilities of responseYN
## 
## pgACC | Predicted |       95% CI
## --------------------------------
##    -6 |      0.29 | [0.11, 0.56]
##    -4 |      0.45 | [0.27, 0.64]
##    -2 |      0.62 | [0.52, 0.70]
##     0 |      0.76 | [0.73, 0.79]
##     2 |      0.86 | [0.80, 0.90]
##     4 |      0.93 | [0.85, 0.96]
##     6 |      0.96 | [0.88, 0.99]
## 
## $vmPFC
## # Predicted probabilities of responseYN
## 
## vmPFC | Predicted |       95% CI
## --------------------------------
##    -4 |      0.89 | [0.77, 0.95]
##    -2 |      0.84 | [0.76, 0.89]
##     0 |      0.77 | [0.74, 0.79]
##     2 |      0.68 | [0.60, 0.76]
##     4 |      0.58 | [0.39, 0.75]
##     6 |      0.47 | [0.22, 0.74]
## 
## $VS
## # Predicted probabilities of responseYN
## 
## VS | Predicted |       95% CI
## -----------------------------
## -4 |      0.73 | [0.62, 0.82]
## -2 |      0.75 | [0.69, 0.80]
##  0 |      0.76 | [0.73, 0.79]
##  2 |      0.77 | [0.72, 0.82]
##  4 |      0.79 | [0.70, 0.86]
##  6 |      0.80 | [0.66, 0.89]
## 
## attr(,"class")
## [1] "ggalleffects" "list"        
## attr(,"model.name")
## [1] "model_2_rt"

pgACC

ggeffect(model_2_rt, terms = c("construct", "pgACC [0, 1]")) %>%
  data.frame() %>%
  mutate(`b [95% CI]` = sprintf("%.2f [%0.2f, %.2f]", predicted, conf.low, conf.high)) %>%
  kable() %>%
  kableExtra::kable_styling()
x predicted std.error conf.low conf.high group b [95% CI]
self-oriented well-being 0.9013345 0.1186334 0.8786410 0.9201699 0 0.90 [0.88, 0.92]
social well-being 0.9055447 0.1222800 0.8829581 0.9241470 0 0.91 [0.88, 0.92]
ill-being 0.2655698 0.0838974 0.2347550 0.2988497 0 0.27 [0.23, 0.30]
self-oriented well-being 0.9346223 0.2305713 0.9009700 0.9573798 1 0.93 [0.90, 0.96]
social well-being 0.9413741 0.2312392 0.9107595 0.9619252 1 0.94 [0.91, 0.96]
ill-being 0.2773118 0.1447915 0.2241538 0.3375907 1 0.28 [0.22, 0.34]

vmPFC

ggeffect(model_2_rt, terms = c("construct", "vmPFC [0, 1]")) %>%
  data.frame() %>%
  mutate(`b [95% CI]` = sprintf("%.2f [%0.2f, %.2f]", predicted, conf.low, conf.high)) %>%
  kable() %>%
  kableExtra::kable_styling()
x predicted std.error conf.low conf.high group b [95% CI]
self-oriented well-being 0.9078067 0.1249982 0.8851522 0.9263642 0 0.91 [0.89, 0.93]
social well-being 0.9146557 0.1344663 0.8917078 0.9331056 0 0.91 [0.89, 0.93]
ill-being 0.2557264 0.0861418 0.2249350 0.2891602 0 0.26 [0.22, 0.29]
self-oriented well-being 0.8691456 0.1913317 0.8203060 0.9062283 1 0.87 [0.82, 0.91]
social well-being 0.8515878 0.1704104 0.8042561 0.8890533 1 0.85 [0.80, 0.89]
ill-being 0.3310591 0.1367623 0.2745920 0.3928496 1 0.33 [0.27, 0.39]

VS

ggeffect(model_2_rt, terms = c("construct", "VS [0, 1]")) %>%
  data.frame() %>%
  mutate(`b [95% CI]` = sprintf("%.2f [%0.2f, %.2f]", predicted, conf.low, conf.high)) %>%
  kable() %>%
  kableExtra::kable_styling()
x predicted std.error conf.low conf.high group b [95% CI]
self-oriented well-being 0.9008684 0.1201313 0.8777676 0.9200013 0 0.90 [0.88, 0.92]
social well-being 0.9059780 0.1228265 0.8833714 0.9245769 0 0.91 [0.88, 0.92]
ill-being 0.2688924 0.0857897 0.2371454 0.3032002 0 0.27 [0.24, 0.30]
self-oriented well-being 0.9102403 0.1535243 0.8824326 0.9319778 1 0.91 [0.88, 0.93]
social well-being 0.9119746 0.1598452 0.8833657 0.9340897 1 0.91 [0.88, 0.93]
ill-being 0.2551545 0.1038266 0.2184360 0.2957100 1 0.26 [0.22, 0.30]

plot predicted effects

by construct

raw_values = data_complete_mod %>%
  select(subjectID, construct, RT, pgACC, vmPFC, VS) %>%
 gather(roi, x, pgACC, VS, vmPFC)

vals = seq(round(min(raw_values$x),1), round(max(raw_values$x), 1), .5)

predicted_sub = ggpredict(model_2_rt, terms = c("pgACC [vals]", "construct", "subjectID"), type = "random") %>%
  data.frame() %>%
  mutate(roi = "pgACC") %>%
  bind_rows(ggpredict(model_2_rt, terms = c("VS [vals]", "construct", "subjectID"), type = "random") %>%
    data.frame() %>%
    mutate(roi = "VS")) %>%
  bind_rows(ggpredict(model_2_rt, terms = c("vmPFC [vals]", "construct", "subjectID"), type = "random") %>%
    data.frame() %>%
    mutate(roi = "vmPFC"))

predicted = ggeffect(model_2_rt, terms = c("pgACC [vals]", "construct")) %>%
  data.frame() %>%
  mutate(roi = "pgACC") %>%
  bind_rows(ggeffect(model_2_rt, terms = c("VS [vals]", "construct")) %>%
    data.frame() %>%
    mutate(roi = "VS")) %>%
  bind_rows(ggeffect(model_2_rt, terms = c("vmPFC [vals]", "construct")) %>%
    data.frame() %>%
    mutate(roi = "vmPFC"))

predicted %>%
  ggplot(aes(x, predicted)) +
  geom_line(data = predicted_sub, aes(color = group, group = interaction(facet, group)), size = .25, alpha = .5) +
  geom_line(aes(color = group), size = 2) +
  facet_grid(~roi) +
  scale_color_manual(values = constructs, name = "") + 
  scale_fill_manual(values = constructs, name = "") + 
  scale_x_continuous(breaks = seq(-6, 6, 2)) +
  labs(x = "\nparameter estimate (SD)", y = "probability of responding yes\n") + 
  plot_aes

by ROI

predicted %>%
  ggplot(aes(x, predicted)) +
  geom_line(data = predicted_sub, aes(color = roi, group = interaction(facet, roi)), size = .25, alpha = .5) +
  geom_line(aes(color = roi), size = 2) +
  facet_grid(~group) +
  scale_color_manual(values = rois, name = "") + 
  scale_fill_manual(values = rois, name = "") + 
  scale_x_continuous(breaks = seq(-6, 6, 2)) +
  labs(x = "\nparameter estimate (SD)", y = "probability of responding yes\n") + 
  plot_aes

RT

vals = seq(round(min(data_complete_mod$RT),1), round(max(data_complete_mod$RT), 1), .1)

predicted_sub = ggpredict(model_2_rt, terms = c("RT [vals]", "construct", "subjectID"), type = "random") %>%
  data.frame()

ggeffect(model_2_rt, terms = c("RT [vals]", "construct")) %>%
  data.frame() %>%
  ggplot(aes(x, predicted)) +
  geom_line(data = predicted_sub, aes(color = group, group = interaction(facet, group)), size = .25, alpha = .5) +
  geom_line(aes(color = group), size = 1.5) +
  scale_color_manual(values = constructs, name = "") + 
  scale_fill_manual(values = constructs, name = "") + 
  labs(x = "\nlog-transformed grand-mean centered RT (s)", y = "probability of responding yes\n") + 
  plot_aes

no RT model

These are the original preregistered analyses not controlling for reaction time, reported in supplementary material.

run model

model_2 = glmer(responseYN ~ construct + pgACC*construct + vmPFC*construct + VS*construct +
                     (1 | subjectID), 
                family = binomial, 
                data = data_complete_mod, control = glmerControl(optimizer = "bobyqa"))

model summary

summary(model_2)
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: binomial  ( logit )
## Formula: responseYN ~ construct + pgACC * construct + vmPFC * construct +  
##     VS * construct + (1 | subjectID)
##    Data: data_complete_mod
## Control: glmerControl(optimizer = "bobyqa")
## 
##      AIC      BIC   logLik deviance df.resid 
##   3355.6   3436.4  -1664.8   3329.6     3681 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.6442 -0.5280  0.3366  0.4098  3.0378 
## 
## Random effects:
##  Groups    Name        Variance Std.Dev.
##  subjectID (Intercept) 0.1913   0.4373  
## Number of obs: 3694, groups:  subjectID, 105
## 
## Fixed effects:
##                                  Estimate Std. Error z value
## (Intercept)                       2.01006    0.10105  19.893
## constructsocial well-being       -0.09374    0.12688  -0.739
## constructill-being               -3.05942    0.11701 -26.146
## pgACC                             0.24530    0.15334   1.600
## vmPFC                            -0.22700    0.15458  -1.469
## VS                                0.12710    0.09061   1.403
## constructsocial well-being:pgACC  0.01747    0.21058   0.083
## constructill-being:pgACC         -0.19387    0.19550  -0.992
## constructsocial well-being:vmPFC -0.07919    0.21652  -0.366
## constructill-being:vmPFC          0.59383    0.19900   2.984
## constructsocial well-being:VS    -0.07809    0.12521  -0.624
## constructill-being:VS            -0.19412    0.12025  -1.614
##                                              Pr(>|z|)    
## (Intercept)                      < 0.0000000000000002 ***
## constructsocial well-being                    0.46005    
## constructill-being               < 0.0000000000000002 ***
## pgACC                                         0.10966    
## vmPFC                                         0.14196    
## VS                                            0.16071    
## constructsocial well-being:pgACC              0.93388    
## constructill-being:pgACC                      0.32137    
## constructsocial well-being:vmPFC              0.71457    
## constructill-being:vmPFC                      0.00285 ** 
## constructsocial well-being:VS                 0.53281    
## constructill-being:VS                         0.10645    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) cnwll- cnstr- pgACC  vmPFC  VS     cw-:AC c-:ACC cw-:PF
## cnstrctwll- -0.627                                                        
## cnstrctll-b -0.725  0.543                                                 
## pgACC        0.210 -0.164 -0.182                                          
## vmPFC       -0.179  0.139  0.155 -0.804                                   
## VS          -0.079  0.060  0.068 -0.087 -0.197                            
## cwll-bn:ACC -0.156  0.261  0.137 -0.718  0.577  0.062                     
## cnstrc-:ACC -0.161  0.129  0.184 -0.773  0.621  0.067  0.562              
## cwll-bn:PFC  0.128 -0.286 -0.112  0.565 -0.703  0.139 -0.799 -0.443       
## cnstrc-:PFC  0.142 -0.109 -0.188  0.615 -0.766  0.151 -0.448 -0.793  0.546
## cwll-bng:VS  0.053 -0.057 -0.045  0.060  0.140 -0.712 -0.075 -0.049 -0.215
## cnstrct-:VS  0.058 -0.044 -0.128  0.066  0.145 -0.745 -0.048 -0.078 -0.104
##             c-:PFC cw-:VS
## cnstrctwll-              
## cnstrctll-b              
## pgACC                    
## vmPFC                    
## VS                       
## cwll-bn:ACC              
## cnstrc-:ACC              
## cwll-bn:PFC              
## cnstrc-:PFC              
## cwll-bng:VS -0.109       
## cnstrct-:VS -0.208  0.537

model table

model_table(model_2)
term b [95% CI] z p
Intercept (self-oriented well-being) 2.01 [1.81, 2.21] 19.89 < .001
Construct (social well-being) -0.09 [-0.34, 0.15] 0.74 .460
Construct (ill-being) -3.06 [-3.29, -2.83] 26.15 < .001
pgACC 0.25 [-0.06, 0.55] 1.60 .110
vmPFC -0.23 [-0.53, 0.08] 1.47 .142
VS 0.13 [-0.05, 0.30] 1.40 .161
Construct (social well-being) x pgACC 0.02 [-0.40, 0.43] 0.08 .934
Construct (ill-being) x pgACC -0.19 [-0.58, 0.19] 0.99 .321
Construct (social well-being) x vmPFC -0.08 [-0.50, 0.35] 0.37 .715
Construct (ill-being) x vmPFC 0.59 [0.20, 0.98] 2.98 .003
Construct (social well-being) x VS -0.08 [-0.32, 0.17] 0.62 .533
Construct (ill-being) x VS -0.19 [-0.43, 0.04] 1.61 .106

marginal probabilities

all

ggeffect(model_2)
## $construct
## # Predicted probabilities of responseYN
## 
## construct                | Predicted |       95% CI
## ---------------------------------------------------
## self-oriented well-being |      0.88 | [0.86, 0.90]
## social well-being        |      0.87 | [0.85, 0.89]
## ill-being                |      0.27 | [0.24, 0.30]
## 
## $pgACC
## # Predicted probabilities of responseYN
## 
## pgACC | Predicted |       95% CI
## --------------------------------
##    -6 |      0.46 | [0.24, 0.70]
##    -4 |      0.55 | [0.39, 0.71]
##    -2 |      0.64 | [0.56, 0.72]
##     0 |      0.72 | [0.70, 0.75]
##     2 |      0.79 | [0.73, 0.84]
##     4 |      0.85 | [0.74, 0.91]
##     6 |      0.89 | [0.75, 0.96]
## 
## $vmPFC
## # Predicted probabilities of responseYN
## 
## vmPFC | Predicted |       95% CI
## --------------------------------
##    -4 |      0.77 | [0.62, 0.87]
##    -2 |      0.75 | [0.67, 0.81]
##     0 |      0.73 | [0.70, 0.75]
##     2 |      0.70 | [0.63, 0.77]
##     4 |      0.68 | [0.53, 0.80]
##     6 |      0.65 | [0.42, 0.84]
## 
## $VS
## # Predicted probabilities of responseYN
## 
## VS | Predicted |       95% CI
## -----------------------------
## -4 |      0.69 | [0.59, 0.78]
## -2 |      0.71 | [0.65, 0.76]
##  0 |      0.72 | [0.70, 0.75]
##  2 |      0.74 | [0.69, 0.78]
##  4 |      0.75 | [0.67, 0.82]
##  6 |      0.76 | [0.64, 0.85]
## 
## attr(,"class")
## [1] "ggalleffects" "list"        
## attr(,"model.name")
## [1] "model_2"

pgACC

ggeffect(model_2, terms = c("construct", "pgACC [0, 1]")) %>%
  data.frame() %>%
  mutate(`b [95% CI]` = sprintf("%.2f [%0.2f, %.2f]", predicted, conf.low, conf.high)) %>%
  kable() %>%
  kableExtra::kable_styling()
x predicted std.error conf.low conf.high group b [95% CI]
self-oriented well-being 0.8811189 0.0989369 0.8592598 0.8999781 0 0.88 [0.86, 0.90]
social well-being 0.8677987 0.0961826 0.8446340 0.8879677 0 0.87 [0.84, 0.89]
ill-being 0.2670746 0.0795909 0.2376673 0.2986948 0 0.27 [0.24, 0.30]
self-oriented well-being 0.9045098 0.1837014 0.8685640 0.9314013 1 0.90 [0.87, 0.93]
social well-being 0.8951448 0.1805836 0.8569846 0.9240255 1 0.90 [0.86, 0.92]
ill-being 0.2772620 0.1410630 0.2253839 0.3359031 1 0.28 [0.23, 0.34]

vmPFC

ggeffect(model_2, terms = c("construct", "vmPFC [0, 1]")) %>%
  data.frame() %>%
  mutate(`b [95% CI]` = sprintf("%.2f [%0.2f, %.2f]", predicted, conf.low, conf.high)) %>%
  kable() %>%
  kableExtra::kable_styling()
x predicted std.error conf.low conf.high group b [95% CI]
self-oriented well-being 0.8855238 0.1026127 0.8635030 0.9043853 0 0.89 [0.86, 0.90]
social well-being 0.8739828 0.1036531 0.8498577 0.8947116 0 0.87 [0.85, 0.89]
ill-being 0.2570822 0.0816078 0.2277362 0.2887955 0 0.26 [0.23, 0.29]
self-oriented well-being 0.8604231 0.1619427 0.8177841 0.8943731 1 0.86 [0.82, 0.89]
social well-being 0.8362314 0.1466111 0.7929990 0.8718932 1 0.84 [0.79, 0.87]
ill-being 0.3330630 0.1345046 0.2772812 0.3939497 1 0.33 [0.28, 0.39]

VS

ggeffect(model_2, terms = c("construct", "VS [0, 1]")) %>%
  data.frame() %>%
  mutate(`b [95% CI]` = sprintf("%.2f [%0.2f, %.2f]", predicted, conf.low, conf.high)) %>%
  kable() %>%
  kableExtra::kable_styling()
x predicted std.error conf.low conf.high group b [95% CI]
self-oriented well-being 0.8793589 0.0998866 0.8570007 0.8986349 0 0.88 [0.86, 0.90]
social well-being 0.8678097 0.0966799 0.8445186 0.8880742 0 0.87 [0.84, 0.89]
ill-being 0.2701836 0.0813785 0.2399069 0.3027594 0 0.27 [0.24, 0.30]
self-oriented well-being 0.8922061 0.1258366 0.8660923 0.9137346 1 0.89 [0.87, 0.91]
social well-being 0.8733311 0.1246071 0.8437650 0.8979786 1 0.87 [0.84, 0.90]
ill-being 0.2571747 0.1004826 0.2213791 0.2965537 1 0.26 [0.22, 0.30]

plot predicted effects

by construct

raw_values = data_complete_mod %>%
  select(subjectID, construct, RT, pgACC, vmPFC, VS) %>%
 gather(roi, x, pgACC, VS, vmPFC)

vals = seq(round(min(raw_values$x),1), round(max(raw_values$x), 1), .5)

predicted_sub = ggpredict(model_2, terms = c("pgACC [vals]", "construct", "subjectID"), type = "random") %>%
  data.frame() %>%
  mutate(roi = "pgACC") %>%
  bind_rows(ggpredict(model_2, terms = c("VS [vals]", "construct", "subjectID"), type = "random") %>%
    data.frame() %>%
    mutate(roi = "VS")) %>%
  bind_rows(ggpredict(model_2, terms = c("vmPFC [vals]", "construct", "subjectID"), type = "random") %>%
    data.frame() %>%
    mutate(roi = "vmPFC"))

predicted = ggeffect(model_2, terms = c("pgACC [vals]", "construct")) %>%
  data.frame() %>%
  mutate(roi = "pgACC") %>%
  bind_rows(ggeffect(model_2, terms = c("VS [vals]", "construct")) %>%
    data.frame() %>%
    mutate(roi = "VS")) %>%
  bind_rows(ggeffect(model_2, terms = c("vmPFC [vals]", "construct")) %>%
    data.frame() %>%
    mutate(roi = "vmPFC"))

predicted %>%
  ggplot(aes(x, predicted)) +
  geom_line(data = predicted_sub, aes(color = group, group = interaction(facet, group)), size = .25, alpha = .5) +
  geom_line(aes(color = group), size = 2) +
  facet_grid(~roi) +
  scale_color_manual(values = constructs, name = "") + 
  scale_fill_manual(values = constructs, name = "") + 
  scale_x_continuous(breaks = seq(-6, 6, 2)) +
  labs(x = "\nparameter estimate (SD)", y = "probability of responding yes\n") + 
  plot_aes

by ROI

predicted %>%
  ggplot(aes(x, predicted)) +
  geom_line(data = predicted_sub, aes(color = roi, group = interaction(facet, roi)), size = .25, alpha = .5) +
  geom_line(aes(color = roi), size = 2) +
  facet_grid(~group) +
  scale_color_manual(values = rois, name = "") + 
  scale_fill_manual(values = rois, name = "") + 
  scale_x_continuous(breaks = seq(-6, 6, 2)) +
  labs(x = "\nparameter estimate (SD)", y = "probability of responding yes\n") + 
  plot_aes

item random effect model

This is a model fit in response to reviewers, included in supplementary material.

run model

model_2_rt_item = glmer(responseYN ~ RT*construct + pgACC*construct + vmPFC*construct + VS*construct +
                     (1 | subjectID) + (1 | item), 
                family = binomial, 
                data = data_complete_mod, control = glmerControl(optimizer = "bobyqa"))

model summary

summary(model_2_rt_item)
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: binomial  ( logit )
## Formula: responseYN ~ RT * construct + pgACC * construct + vmPFC * construct +  
##     VS * construct + (1 | subjectID) + (1 | item)
##    Data: data_complete_mod
## Control: glmerControl(optimizer = "bobyqa")
## 
##      AIC      BIC   logLik deviance df.resid 
##   2786.0   2891.6  -1376.0   2752.0     3677 
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -10.9511  -0.3244   0.1677   0.3614   4.6525 
## 
## Random effects:
##  Groups    Name        Variance Std.Dev.
##  subjectID (Intercept) 0.3501   0.5917  
##  item      (Intercept) 1.0609   1.0300  
## Number of obs: 3694, groups:  subjectID, 105; item, 48
## 
## Fixed effects:
##                                  Estimate Std. Error z value
## (Intercept)                       1.96805    0.30826   6.384
## RT                               -3.51299    0.31973 -10.988
## constructsocial well-being        0.41306    0.43415   0.951
## constructill-being               -3.24427    0.41970  -7.730
## pgACC                             0.39922    0.19094   2.091
## vmPFC                            -0.36949    0.19099  -1.935
## VS                                0.13200    0.11471   1.151
## RT:constructsocial well-being     0.34789    0.42647   0.816
## RT:constructill-being             4.05498    0.39109  10.368
## constructsocial well-being:pgACC  0.13295    0.26468   0.502
## constructill-being:pgACC         -0.41633    0.23229  -1.792
## constructsocial well-being:vmPFC -0.22625    0.27384  -0.826
## constructill-being:vmPFC          0.66764    0.23558   2.834
## constructsocial well-being:VS    -0.07297    0.16141  -0.452
## constructill-being:VS            -0.12616    0.14516  -0.869
##                                              Pr(>|z|)    
## (Intercept)                        0.0000000001721625 ***
## RT                               < 0.0000000000000002 ***
## constructsocial well-being                     0.3414    
## constructill-being                 0.0000000000000107 ***
## pgACC                                          0.0365 *  
## vmPFC                                          0.0530 .  
## VS                                             0.2498    
## RT:constructsocial well-being                  0.4146    
## RT:constructill-being            < 0.0000000000000002 ***
## constructsocial well-being:pgACC               0.6155    
## constructill-being:pgACC                       0.0731 .  
## constructsocial well-being:vmPFC               0.4087    
## constructill-being:vmPFC                       0.0046 ** 
## constructsocial well-being:VS                  0.6512    
## constructill-being:VS                          0.3848    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation matrix not shown by default, as p = 15 > 12.
## Use print(x, correlation=TRUE)  or
##     vcov(x)        if you need it

model table

model_table(model_2_rt_item)
term b [95% CI] z p
Intercept (self-oriented well-being) 1.97 [1.36, 2.57] 6.38 < .001
RT -3.51 [-4.14, -2.89] 10.99 < .001
Construct (social well-being) 0.41 [-0.44, 1.26] 0.95 .341
Construct (ill-being) -3.24 [-4.07, -2.42] 7.73 < .001
pgACC 0.40 [0.02, 0.77] 2.09 .037
vmPFC -0.37 [-0.74, 0.00] 1.93 .053
VS 0.13 [-0.09, 0.36] 1.15 .250
RT x Construct (social well-being) 0.35 [-0.49, 1.18] 0.82 .415
RT x Construct (ill-being) 4.05 [3.29, 4.82] 10.37 < .001
Construct (social well-being) x pgACC 0.13 [-0.39, 0.65] 0.50 .615
Construct (ill-being) x pgACC -0.42 [-0.87, 0.04] 1.79 .073
Construct (social well-being) x vmPFC -0.23 [-0.76, 0.31] 0.83 .409
Construct (ill-being) x vmPFC 0.67 [0.21, 1.13] 2.83 .005
Construct (social well-being) x VS -0.07 [-0.39, 0.24] 0.45 .651
Construct (ill-being) x VS -0.13 [-0.41, 0.16] 0.87 .385

marginal probabilities

all

ggeffect(model_2_rt_item)
## $RT
## # Predicted probabilities of responseYN
## 
##    RT | Predicted |       95% CI
## --------------------------------
## -1.50 |      0.98 | [0.97, 0.99]
## -1.00 |      0.96 | [0.93, 0.97]
## -0.50 |      0.88 | [0.84, 0.92]
##  0.00 |      0.73 | [0.66, 0.80]
##  0.50 |      0.50 | [0.40, 0.60]
##  1.00 |      0.26 | [0.18, 0.37]
## 
## $construct
## # Predicted probabilities of responseYN
## 
## construct                | Predicted |       95% CI
## ---------------------------------------------------
## self-oriented well-being |      0.91 | [0.85, 0.95]
## social well-being        |      0.94 | [0.89, 0.96]
## ill-being                |      0.22 | [0.13, 0.33]
## 
## $pgACC
## # Predicted probabilities of responseYN
## 
## pgACC | Predicted |       95% CI
## --------------------------------
##    -6 |      0.36 | [0.14, 0.66]
##    -4 |      0.50 | [0.30, 0.71]
##    -2 |      0.65 | [0.52, 0.76]
##     0 |      0.77 | [0.70, 0.83]
##     2 |      0.86 | [0.79, 0.92]
##     4 |      0.92 | [0.83, 0.97]
##     6 |      0.96 | [0.86, 0.99]
## 
## $vmPFC
## # Predicted probabilities of responseYN
## 
## vmPFC | Predicted |       95% CI
## --------------------------------
##    -4 |      0.90 | [0.77, 0.96]
##    -2 |      0.85 | [0.76, 0.91]
##     0 |      0.78 | [0.71, 0.84]
##     2 |      0.70 | [0.58, 0.79]
##     4 |      0.60 | [0.38, 0.78]
##     6 |      0.49 | [0.21, 0.77]
## 
## $VS
## # Predicted probabilities of responseYN
## 
## VS | Predicted |       95% CI
## -----------------------------
## -4 |      0.73 | [0.58, 0.83]
## -2 |      0.75 | [0.66, 0.83]
##  0 |      0.77 | [0.70, 0.83]
##  2 |      0.80 | [0.72, 0.86]
##  4 |      0.82 | [0.71, 0.89]
##  6 |      0.84 | [0.69, 0.92]
## 
## attr(,"class")
## [1] "ggalleffects" "list"        
## attr(,"model.name")
## [1] "model_2_rt_item"

pgACC

ggeffect(model_2_rt_item, terms = c("construct", "pgACC [0, 1]")) %>%
  data.frame() %>%
  mutate(`b [95% CI]` = sprintf("%.2f [%0.2f, %.2f]", predicted, conf.low, conf.high)) %>%
  kable() %>%
  kableExtra::kable_styling()
x predicted std.error conf.low conf.high group b [95% CI]
self-oriented well-being 0.9116221 0.3103964 0.8488022 0.9498825 0 0.91 [0.85, 0.95]
social well-being 0.9346029 0.3199374 0.8841754 0.9639701 0 0.93 [0.88, 0.96]
ill-being 0.2155277 0.2958810 0.1333292 0.3291558 0 0.22 [0.13, 0.33]
self-oriented well-being 0.9389360 0.3698387 0.8816297 0.9694598 1 0.94 [0.88, 0.97]
social well-being 0.9605250 0.3804260 0.9202824 0.9808747 1 0.96 [0.92, 0.98]
ill-being 0.2126499 0.3240925 0.1251833 0.3376428 1 0.21 [0.13, 0.34]

vmPFC

ggeffect(model_2_rt_item, terms = c("construct", "vmPFC [0, 1]")) %>%
  data.frame() %>%
  mutate(`b [95% CI]` = sprintf("%.2f [%0.2f, %.2f]", predicted, conf.low, conf.high)) %>%
  kable() %>%
  kableExtra::kable_styling()
x predicted std.error conf.low conf.high group b [95% CI]
self-oriented well-being 0.9170690 0.3129186 0.8569107 0.9533132 0 0.92 [0.86, 0.95]
social well-being 0.9408932 0.3251703 0.8938003 0.9678544 0 0.94 [0.89, 0.97]
ill-being 0.2081341 0.2967131 0.1281119 0.3198061 0 0.21 [0.13, 0.32]
self-oriented well-being 0.8842881 0.3460413 0.7950193 0.9377256 1 0.88 [0.80, 0.94]
social well-being 0.8976827 0.3455053 0.8167617 0.9452626 1 0.90 [0.82, 0.95]
ill-being 0.2615243 0.3191904 0.1592709 0.3983226 1 0.26 [0.16, 0.40]

VS

ggeffect(model_2_rt_item, terms = c("construct", "VS [0, 1]")) %>%
  data.frame() %>%
  mutate(`b [95% CI]` = sprintf("%.2f [%0.2f, %.2f]", predicted, conf.low, conf.high)) %>%
  kable() %>%
  kableExtra::kable_styling()
x predicted std.error conf.low conf.high group b [95% CI]
self-oriented well-being 0.9106727 0.3113140 0.8470583 0.9494078 0 0.91 [0.85, 0.95]
social well-being 0.9351202 0.3203609 0.8849579 0.9642925 0 0.94 [0.88, 0.96]
ill-being 0.2152117 0.2968624 0.1328915 0.3291677 0 0.22 [0.13, 0.33]
self-oriented well-being 0.9208444 0.3249766 0.8601977 0.9565124 1 0.92 [0.86, 0.96]
social well-being 0.9386109 0.3371729 0.8875865 0.9673279 1 0.94 [0.89, 0.97]
ill-being 0.2161998 0.3029340 0.1321948 0.3330969 1 0.22 [0.13, 0.33]

plot predicted effects

by construct

raw_values = data_complete_mod %>%
  select(subjectID, construct, RT, pgACC, vmPFC, VS) %>%
 gather(roi, x, pgACC, VS, vmPFC)

vals = seq(round(min(raw_values$x),1), round(max(raw_values$x), 1), .5)

predicted_sub = ggpredict(model_2_rt_item, terms = c("pgACC [vals]", "construct", "subjectID"), type = "random") %>%
  data.frame() %>%
  mutate(roi = "pgACC") %>%
  bind_rows(ggpredict(model_2_rt_item, terms = c("VS [vals]", "construct", "subjectID"), type = "random") %>%
    data.frame() %>%
    mutate(roi = "VS")) %>%
  bind_rows(ggpredict(model_2_rt_item, terms = c("vmPFC [vals]", "construct", "subjectID"), type = "random") %>%
    data.frame() %>%
    mutate(roi = "vmPFC"))

predicted = ggeffect(model_2_rt_item, terms = c("pgACC [vals]", "construct")) %>%
  data.frame() %>%
  mutate(roi = "pgACC") %>%
  bind_rows(ggeffect(model_2_rt_item, terms = c("VS [vals]", "construct")) %>%
    data.frame() %>%
    mutate(roi = "VS")) %>%
  bind_rows(ggeffect(model_2_rt_item, terms = c("vmPFC [vals]", "construct")) %>%
    data.frame() %>%
    mutate(roi = "vmPFC"))

predicted %>%
  ggplot(aes(x, predicted)) +
  geom_line(data = predicted_sub, aes(color = group, group = interaction(facet, group)), size = .25, alpha = .5) +
  geom_line(aes(color = group), size = 2) +
  facet_grid(~roi) +
  scale_color_manual(values = constructs, name = "") + 
  scale_fill_manual(values = constructs, name = "") + 
  scale_x_continuous(breaks = seq(-6, 6, 2)) +
  labs(x = "\nparameter estimate (SD)", y = "probability of responding yes\n") + 
  plot_aes

by ROI

predicted %>%
  ggplot(aes(x, predicted)) +
  geom_line(data = predicted_sub, aes(color = roi, group = interaction(facet, roi)), size = .25, alpha = .5) +
  geom_line(aes(color = roi), size = 2) +
  facet_grid(~group) +
  scale_color_manual(values = rois, name = "") + 
  scale_fill_manual(values = rois, name = "") + 
  scale_x_continuous(breaks = seq(-6, 6, 2)) +
  labs(x = "\nparameter estimate (SD)", y = "probability of responding yes\n") + 
  plot_aes

RT

vals = seq(round(min(data_complete_mod$RT),1), round(max(data_complete_mod$RT), 1), .1)

predicted_sub = ggpredict(model_2_rt_item, terms = c("RT [vals]", "construct", "subjectID"), type = "random") %>%
  data.frame()

ggeffect(model_2_rt_item, terms = c("RT [vals]", "construct")) %>%
  data.frame() %>%
  ggplot(aes(x, predicted)) +
  geom_line(data = predicted_sub, aes(color = group, group = interaction(facet, group)), size = .25, alpha = .5) +
  geom_line(aes(color = group), size = 1.5) +
  scale_color_manual(values = constructs, name = "") + 
  scale_fill_manual(values = constructs, name = "") + 
  labs(x = "\nlog-transformed grand-mean centered RT (s)", y = "probability of responding yes\n") + 
  plot_aes